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Streamers are a mode of dielectric breakdown of a gas in a strong electric field: A sharp nonlinear 
ionization wave propagates into a non-ionized gas, leaving a nonequilibrium plasma behind. The 
ionization avalanche in the tip of the wave is due to free electrons being accelerated in the strong field 
and ionizing the gas by impact. This chain reaction deeper in the wave is suppressed by the generated 
free charges screening the field. Simulations of streamers show two widely separated spatial scales: 
the width of the charged layer where the electron density gradients and the ionization rate are very 
large (O(fim)), and the width of the electrically screened, finger-shaped and ionized region (O(mm)). 
We thus recently have suggested to analyze first the properties of the charge-ionization-layer on the 
inner scale on which it is almost planar, and then to understand the streamer shape on the outer 
scale as the motion of an effective interface as is done in other examples of nonequilibrium pattern 
formation. The first step thus is the analysis of the inner dynamics of planar streamer fronts. For 
these, we resolve the long-standing question about what determines the front speed, by applying 
the modern insights of pattern formation to the streamer equations used in the recent simulations. 
These include field-driven impact ionization, electron drift and diffusion and the Poisson equation 
for the electric field. First, in appropriately chosen dimensionless units only one parameter remains 
to characterize the gas, the dimensionless electron diffusion constant D; for typical gases under 
normal conditions D « 0.1-0.3. Then we determine essentially all relevant properties of planar 
streamer fronts. Technically, we identify the propagation of streamer fronts as an example of front 
propagation into unstable states. In terms of the marginal stability scenario we then find: the front 
approached asymptotically starting from any sufficiently localized initial condition (the "selected 
front"), is the steepest uniformly translating front solution, which is physical and stable. Negatively 
charged fronts are selected by linear marginal stability, which allows us to derive their velocity 
analytically. Positively charged fronts can only propagate due to electron diffusion against the 
electric field; as a result their behavior is singular in the limit of D — > 0. For D < 1, these fronts 
are selected by nonlinear marginal stability and we have to apply numerical methods for predicting 
the selected front velocity. For larger D, linear marginal stability applies and the velocity can be 
determined analytically. Numerical integrations of the temporal evolution of planar fronts out of 
localized initial conditions confirm all our analytical and numerical predictions for the selection. 
Finally, our general predictions for the selected front velocity and for the degree of ionization of 
the plasma are in semi-quantitative agreement with recent numerical solutions of three-dimensional 
streamer propagation. This gives credence to our suggestion, that the front analysis on the inner 
()im) scale yields the moving boundary conditions for a moving "streamer interface" , whose pattern 
formation is governed by the evolution of the fields on the outer (mm) scale. 
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I. INTRODUCTION 

Discharges are nonequilibrium ionization processes oc- 
curring in initially non-ionized matter exposed to strong 
electric fields. Depending on the spatio-temporal char- 
acteristics of the electric field and on the ionization and 
charge transport properties of the medium, discharges 
can assume many different modes of appearance. In par- 
ticular in gases under approximately normal conditions 
one distinguishes phenomenologically between stationary 
modes like arc, glow or dark discharges and transient 
phenomena like leaders, the initial stages of sparks, and 
streamers jj]-[§, which occur e.g. in silent discharges 0. 
The latter nonstationary discharges often form the initial 
state of a discharge that later on becomes stationary. We 



will focus here on an essential element of many transient 
discharge phenomena, the initial field-driven ionization 
wave. 

The conceptually simplest problem of this kind has be- 
come known as the streamer problem in a non-attaching 
gas. It treats the dynamics of the free electrons and pos- 
itive ions in a homogeneous gas at rest taking the follow- 
ing mechanisms into account: (i) impact ionization, the 
process in which a free electron accelerated in a strong 
local field ionizes a neutral molecule, generating a new 
free electron and a positive ion; (ii) drift and diffusion 
of charged particles, in particular of the electrons whose 
mobility is much larger than that of the ions; (Hi) the 
coupling of the electric field to the charges through the 
Poisson equation of electrostatics. 



Recent numerical simulations [|3|J9[] of a basic model in- 
corporating these physical ingredients for parameter val- 
ues appropriate for nitrogen under normal conditions re- 
veal that a streamer consists of a sharp nonlinear ion- 
ization front which propagates into a non-ionized gas, 
leaving a weakly ionized nonequilibrium plasma behind. 
The underlying mechanism is that in the leading edge 
of the front the electrons are accelerated by the large 
imposed electric field; this causes the build-up of an elec- 
tron avalanche due to impact ionization. The generated 
free charges eventually screen the field and thus suppress 
further ionization. It is the nonlinear balance between 
these two nonequilibrium processes, namely the ioniza- 
tion avalanche and the electric screening, which deter- 
mines the dynamics of the ionization front and the state 
of the plasma behind it. In confined geometries, stream- 
ers usually have a nontrivial fingerlike shape, as is illus- 
trated by the snapshots in Fig. 1 of streamer dynamics 
taken from the simulations of Vitello et al. |9). As the 
sharpness of the electron density profiles in Fig. 1 illus- 
trates, the "passive body" of the finger is separated from 
the external non-ionized gas by a very narrow region — 
of width of order microns — in which essentially all the 
action is occurring. This width has to be compared to 
the size of the filament, which is of order millimeters. It 
is in this narrow layer that most of the ionization process 
is taking place. In this same region, there is a nonzero 
charge density, and consequently, also a very large elec- 
tric field gradient. These features indicate that there are 
two different spatial scales in this process, an "inner" 
scale associated with the thickness of the zone where the 
ionization takes place, and an "outer" one where the spa- 
tial variations are set by the size of the finger and the 
external experimental geometry. It is precisely for these 
reasons that accurate simulations are extremely demand- 
ing and that they were accomplished only recently by 
Dhali and Williams || and by Vitello et al. [|| (See also 

Such a separation of scales is strongly reminiscent of 
what occurs in combustion fronts jO|L2 ]. A combus- 
tion front is a narrow layer of thickness £i n to which 
the combustion is essentially confined, while outside of 
it, the temperature field varies on a much longer scale 
lout- Physically, such sharp combustion fronts occur in 
the limit when the chemical reaction rates involved in 
the combustion are very fast once a sufficient tempera- 
ture is reached. It has been shown that, on the basis of an 
asymptotic expansion to lowest order in the small param- 
eter e = tin I (-out using matched asymptotic expansions 
|L3| , [14| , the problem can be analyzed in terms of the prop- 
agation of an "effective interface" . More specifically, one 
first solves the so-called inner problem of a locally almost 
planar reaction zone. This permits to relate the temper- 
ature and chemical composition fields on both sides of 
the front (at distances L such that < m « L « tout) and 
to determine the local front velocity as a function of local 
curvature and fields. On the scale of the remaining outer 
problem, these relations then play the role of boundary 



conditions and of a kinetic equation for the effective mov- 
ing interface of zero thickness. Besides in combustion, the 
technique of asymptotic matching to obtain an effective 
interface description has also been applied to chemical 
waves pa] , thermal plumes [ fL6[ and to phase field mod- 
els of solidification |L7||J). 

In spite of some important differences between combus- 
tion and streamer fronts as discussed in the appendix, a 
similar approach appears possible for streamers. As dis- 
cussed also in Q, building on such a reduced descrip- 
tion of streamer dynamics appears very desirable, not 
only because it might make numerical studies much eas- 
ier, but also because it will allow us to draw upon the 
knowledge and methods which have been developed in 
the last decade in the field of interfacial pattern forma- 
tion and dynamics |2Q] . The first step towards this goal 
is to determine the field dependence of the velocity and 
the ionization and charge profile of a planar front which 
propagates into the non-ionized region. We thus analyze 
in this paper the inner problem for a planar streamer 
front. This allows us to reduce the problem to effectively 
one dimension. Our analysis clearly identifies the prob- 
lem of streamer front propagation as an example of front 
propagation into unstable states. Physically, the instabil- 
ity of the non-ionized gas against charge fluctuations can 
be traced back to the fact that any small electron density 
gets amplified by the impact ionization. As is standard 
for front propagation into unstable states 1 21 - 25 1 , we find 
that the one-dimensional streamer equations exhibit a 
one-parameter family of uniformly translating front so- 
lutions, parametrized by their velocity. As usual [pl|-p5[, 
the question is then to decide which of these front so- 
lutions is the dynamically selected one, i.e., is the one 
reached at long times after a localized ionized region has 
been created by some initial ionization event. The exist- 
ing knowledge of front propagation into unstable states 
p2| , P3| provides us with an educated guess for the se- 
lected velocity, which we confirm with the help of nu- 
merical studies. Taken together, our results provide an 
essentially complete solution of the inner problem of pla- 
nar streamer fronts. 

In itself, the idea to analyze the planar fronts of a 
streamer model is not new — we refer to p6|-[29| for ear- 
lier work. Apart from the fact, that the authors from the 
seventies [p6|-p8| investigate different models, which are 
more inspired by equilibrium concepts (e.g., the ioniza- 
tion behind the front is determined by thermal ionization, 
where the electron temperature is raised by application 
of strong electric fields), our work casts new light on this 
old problem from two different angles: 
First of all, it was empirically noted, that the standard 
approach to analyze uniformly translating fronts, failed 
to determine a unique propagation velocity, given the 
field and the gas parameters. Turcotte and Ong |2(| 
clearly state this failure of their theory (this "great de- 
fect" of their theory is recalled in Fowler's reviews [p8[ ) 
and suggest, that a unique solution might be determined 
by a dynamical stability analysis. Albright and Tidman 



p7j then perform such a stability analysis, but not in 
a systematic way, and they draw incorrect conclusions. 
D'yakonov and Kachorovskii p9] also find the indetermi- 
nacy of the speed of uniformly translating planar fronts, 
now for an approximated version of our model, and pro- 
pose to solve this by using the tip radius of the streamer 
finger as an extra length scale, which they, however, can- 
not determine. We, in contrast, trace the indeterminacy 
of the velocity from the analysis of uniformly translating 
streamer front solutions to the fact, that this is an exam- 
ple of front propagation into unstable states. Applying 
the concepts explained above, we solve the selection prob- 
lem for planar fronts without additional assumptions or 
approximations. We argue that a particular front solu- 
tion out of a whole family of dynamically stable solutions 
is selected, because it is the only one compatible with the 
initial condition of a localized ionization seed. 
Secondly, this result is the first ingredient for studying the 
formation of patterns, in particular of the tip radius — 
we do not attempt to model global features of the pattern 
formation with our planar front analysis. Our approach 
thus is very different in spirit to the earlier investiga- 
tions: As also stressed in |19[, in an effective interface 
description based on a matched asymptotic expansion, 
the results of weakly curved, almost planar fronts are es- 
sentially used locally everywhere in the interface region: 
They enter the analysis on the outer scale as boundary 
conditions at the moving interface. It is on this outer 
scale that pattern formation problems like the size, veloc- 
ity and shape of the streamer should be analyzed. Once 
our results on planar fronts will be extended to weakly 
curved fronts, all the necessary ingredients to tackle these 
questions appear to be available. 

The main results of our present analysis of the streamer 
equations used in the simulations Mp) can be summa- 
rized as follows: 

(a) Dimensional analysis shows that in dimensionless 
units, a single parameter remains to characterize the gas, 
the dimensionless electron diffu sion coefficient D charac- 
teristic of the gas [see Eq. ( 2.10J )] . For gases under normal 
conditions, D is small, of order 0.1-0.3. 

(b) The length scale set by the electron i mpa ct ion- 
ization coefficient [the coefficient a$ 1 in Eq. ( |2.6| )1 is on 
the order of microns for nitrogen. For D < 1 the thick- 
ness (.in of the charged layer is on the order of this same 
ionization length for negatively charged streamer fronts 
(NSF) Q . Given that typical streamer diameters found 
in the simulations are of the order of 1mm, e = £i n /(out is 
at most of order 10~ 2 ; this justifies an effective interface 
description of streamer dynamics. 

(c) We find that electron diffusion acts as a singu- 
lar perturbation for positively charged streamer fronts 
(PSF): without diffusion, such fronts can not propagate, 
but with any nonzero D, they do. As a result, the be- 
havior is singular in the limit D — > 0: for D — 0(1), the 
thickness £i n is again of order of the ionization length, 
but for D — ► the electron density and its gradients di- 
verge due to the appearance of another smaller length 



scale (of order D/(xq). 

(d) The electron density generated by the propagating 
front is again basically set by dimensional analysis for 
NSF. We calculate for D < 1.5 the dependence of the 
dimensionless electron density er~ behind the front on 
the electric field E + far ahead of our planar front. Our 
results compare favorably with those extracted from the 
simulations |9| , according to the prescriptions of the the- 
ory of matched asymptotic expansions J13LH ■ Namely, 
E + is not the field value at the electrode position, but the 
value obtained by extrapolating the slowly varying outer 
field to the front position. We also calculate the full D 
and E + dependence of the electron density o~~ behind 
the front of PSF for D < 1.5. 

(e) The dynamically relevant ("selected") front veloc- 
ity Vf is a unique function of E + and D. The analysis 
confirms the strong asymmetry between NSF and PSF 
also found in the simulations MM for fronts propagating 
into an essentially non-ionized region. The asymmetry is 
the stronger the smaller D and disappears for D 3> 1. 

(f) For NSF, Vf is given by the so-called linear 
marginal stability velocity v* p2| — see Eq. (5.3) below. 



For parameter values used in the simulations, we find 
that Vf is typically 30 to 40 % higher than the electron 
drift velocity just in front of the streamer head, which 
agrees semi-quantitatively with the findings of Vitello et 
at. |. 

(g) We find that PSF propagate for any nonzero value 
of the dimensionless electron diffusion coefficient D. Due 
to the singular behavior as D — > 0, we find that fronts 
propagate with a unique velocity v^ predicted by the 
so-called nonlinear marginal stability mechanism |23J for 
small D. For the Townsend expression used in the sim- 
ulations HH, this happens below a well-defined field- 
dependent value of D of order unity [see Fig. 3]. Above 
this threshold value, PSF propagate with the linear 
marginal stability value v*. 

In this paper, our main focus will be on those results 
that are of greatest interest from the point of view of 
understanding the generation of low temperature plas- 
mas by the streamer mechanism. We note, however, that 
t he eq uations for planar streamer fronts [Eqs. (3.11) and 
( |3.12|) below] appear to be of interest in their own right. 
As will be discussed briefly in Section V, our streamers 
have several features in common with the celebrated non- 



linear diffusion equation studied in mathematics 
since the early work of Kolmogorov et al. |3^] and Fisher 
fl34f ; at the same time, however, they are sufficiently more 
complicated that they appear to present new challenges 
from a mathematical point of view. 

This paper is organized as follows. In Section II we in- 
troduce the basic equations for streamer formation, and 
perform a dimensional analysis for the inner problem of 
streamer fronts. In Section III, we discuss the stability 
of the basic homogeneous states of interest, the homo- 
geneous non-ionized state and the homogeneous weakly 
ionized state. We also discuss the physical mechanism 
of streamer formation and the proper initial and bound- 



ary conditions to study these in the case of planar fronts, 
which allow us to simplify the equations describing planar 
front dynamics. In Section IV we demonstrate that there 
exists a one-parameter family of uniformly translating 
fronts characterized by a continuous range of front veloc- 
ities v. We also briefly show how in the case D = 0, the 
equations for uniformly translating fronts can be solved 
analytically. These solutions, which turn out to be useful 
as a small-D approximation for NSF, yield an explicit 
expression for the electron density o~ behind the NSF 
in terms of the field E + just ahead of it. This is followed 
by an analysis of the general case D^O; then the equa- 
tions can not be solved analytically, but we demonstrate 
that there still is a one-parameter family of uniformly 
translating front solutions. For PSF, we show that the 
limit D — ► is singular; we discuss this limit in detail 
and show that it accounts for the strong asymmetry be- 
tween PSF and NSF for realistic values oi D. In Section 
V we then summarize some of the main results 21 -2q] 
concerning the so-called selection problem, the question 
which particular front solution from the family is reached 
asymptotically for large times for a large class of initial 
conditions. Application of these concepts allows us to 
predict the shape and velocity of the dynamically rele- 
vant front solution (the selected front) and the value of 
the electron density generated behind it. This yields the 
various selection results for NSF and for PSF, summa- 
rized in points (c)-(g) above, and leads us to predict that 
the behavior of PSF in the limit D — > is singular. In 
Section VI we present numerical simulations of the full 
partial differential equations for planar streamer dynam- 
ics; starting from various initial conditions, we illustrate, 
that in all cases we have studied the long time dynam- 
ics of the system is characterized by a NSF and a PSF 
whose behavior is in full agreement with our predictions. 
In the concluding section we finally reflect on our results 
and on the future steps to be taken to arrive at an ef- 
fective interface description of streamer dynamics. In an 
appendix we discuss differences and similarities between 
combustion and streamer fronts. 



II. MODELING AND DIMENSIONAL ANALYSIS 



A. The Minimal Streamer Model 



For simulating the dynamical development of stream- 
ers out of a macroscopic initial ionization seed in a so- 
called non-attaching gas like N2 under normal conditions, 
Dhali and Williams J8| and Vitello et al. || use the fol- 
lowing set of deterministic continuum equations for the 
electron density n e , the ion density n + and the electric 
field £: balance equations for electrons and ions, 



where the fact that the two source terms are the same 
is due to charge conservation in an ionization event; the 
Poisson equation, 



Vr • £ — — (n+ - Tie) , 

£0 



(2.3) 



and the approximate phcnomcnological expressions 

j e = -n e He £ - D e V R n e , (2.4) 

J+=0, (2.5) 

source = \n e ^ e £\ a e~ E ol\ E \ . (2.6) 

Apart from the fa ct th at we will allow for a slight gen- 
eralization of Eq. (2.6), these are the equations that we 
will investigate analytically below. 

In these equations, j e and j+ are the particle current 
densities of electrons and positive ions, and e is the ab- 
solute value of the electron charge. The (dimensional) 
spatial coordinates are denoted by R, and Vr is the gra- 
dient with respect to these coordinates. The use of only 



Poisson's law of electrostatics, Eq. (2.3), means that all 



d t n e + Vr • j e = source 
d t n+ + Vr • j+ = source 



(2.1) 
(2.2) 



magnetic fields as well as terms in the Maxwell equa- 
tions associated with time-dependences of the fields, are 
neglected 35. 

The ele ctro n particle current density j e is approxi- 
mated in (2.4) as the sum of a drift and a diffusion term. 
Note that this diffusion approximation implies that the 
electron mean free path must be small with respect to 
the scale of variation £j„ of the electric field. This condi- 
tion is just about satisfied for the parameter values taken 
for N2 in the simulations, except possibly at the highest 
field values (see also the discussion in the concluding Sec- 
tion VII). The electron drift velocity is taken to be lin- 
ear in the field £, with /j, e the (positive) electron mobil- 
ity. The electron diffusion coefficient D e and the mobil- 
ity [i e are treated here as independent coefficients, since 
they effectively depend on the field strength || (only in 
the low-field limit are they related by the Einstein re- 
lation). More generally, the diffusion coefficient should 
be replaced by a diffusion tensor, which is diagonal in a 
reference frame with one axis along the electric field. Its 
longitudinal component, the only relevant one for planar 
fronts perpendicular to E, is somewhat smaller than the 
transverse one. Since we will see, that N2 reaches a typi- 
cal degree of ionization of only 10 -5 , density fluctuations 
of the non-ionized gas can be neglected and the mean 
free path of the electrons and therefore /i e and D e can 
be taken as independent of the degree of ionization. 

The ionic current is neglected according to Eq. (|2.5|), 
since the mobility of ions is at least two orders of magni- 
tude smaller than that of the electrons [||. In particular, 
for the analysis of the inner scale, that we will perform 
in the present paper, j + is negligible. 

The source (2.6) finally accounts for the creation of 
free charges by impact ionization. If the product of elec- 
tric field £ and electronic mean free path £ m f p is large 
enough, free electrons can gain sufficient kinetic energy to 



ionize neutral molecules. Accordingly there is a thresh- 
old field \£\ = E Q oc Imfp^ 1 - For \£\ > E Q the proba- 
bility that a scattering event carries at least the ioniza- 
tion energy is large. The effective ionization cross-section 
ct cs (|£|) then essentially saturates, while for \£\ <C £b 
the ionization rate per scattering event is largely sup- 
pressed. The source-term is given by the ionization rate, 
which can be calculated as the product of the drift cur- 
rent of free electrons \n e n e £\ times the target particle 
density n n of the neutral gas times the effective ioniza- 
tion cross-section <t cs (|£|). Commonly, a phenomcnolog- 
ical ionization coefficient a(|£|) = n n a cs (\£\) is used, 
(which clearly has dimension of inverse length,) whose 
field threshold behavior in the Townsend approxim atio n 
a(|£|) = aoexp(— Eq/\£J) |J is expressed by Eq. (2.6). 
As discussed by Raizer H , in the approximation that ev- 
ery collision is ionizing, if the electron carries an energy 
larger than the ionization energy I, we have 

and E w I/(e £ m f P ) ■ 



a 



r 1 

mfp 



(2.7) 

Since in much of our analysis the specific form of o/(|£|) 
is not needed, we w ill use a slightly more general formu- 
lation in Eq. (2.11) below. 

In the source-term, ionization due to the photons 
also created in recombination or scattering events, is 
neglected. This is motivated by the ionization cross- 
sections due to photons being much smaller than those 
due to electrons. Note that, if photo-ionization is taken 
into account, the dynamical equations become nonlocal. 

No sink term needs to be included for the analysis of 
the inner problem, since the recombination length at a 
degree of ionization of order 10 -5 that we will derive 
below is very large as compared with the front width 
li n . (For this reason, the inner problem is the same for 
streamers and leaders H: the difference between these 
discharge modes, which consists in the fact that recom- 
bination is non-negligible in the plasma body of leaders, 
would come into play only when solving, at a later stage, 
the outer problem.) The fact that the degree of ion- 
ization remains small is also the reason that saturation 



effects are neglected in ( |2.6| ). 

In contrast to the situation in N2, that is described 
by our model equations, in attaching gases like O2, a 
third charged species plays a role, namely negative ions 
formed by a neutral molecule catching a free electron. 
For a description of the physics of such attaching gases 
and simulations thereof, see, e.g., J36[ |. 

The equations above are deterministic. Thermal fluc- 
tuations in fact can be neglected, since even an unphys- 
ically small ionization energy of 3 eV leads to a Boltz- 
mann factor of 10~ 52 at room temperature. Also other 
stochastic effects are not accounted for in the simulations 
we compare to. We further discuss possible stochastic ef- 
fects in the experiments in the conclus ion . 

Finally, the dynamical system ( 2. 1 )-(2.£ ) must be com- 
plemented by: 

(i) boundary conditions: as will be discussed in detail 
in Section III, for the problem of front propagation, these 



are specified by the value E + of the electric field far ahead 
of the front, where the total charge density vanishes. 

(ii) initial conditions: we ignore the details of the 
plasma nucleation event (e.g. triggering by radiation 
from an external source), and assume that at t = a 
small well- localized ionization seed is present. The pre- 
cise meaning, for our problem, of "well-localized" will be 
made clear in Section V. 



B. Dimensional Analysis 

In order to identify the physical scales and intrinsic 
parameters of our problem, we reduce Eqs. (2.1)- (2. 6) to 
a dimcnsionless form. The most natural scale of length 
and electric field are the ionization length li on — ctg 
and the threshold field Eq of the ionization rate ([2.6]). 
The velocity scale is then the electron drift velocity at 
this field strength, vo = ^ e Eo, leading to a time unit 
to = (ctoMe-E'o) -1 , an d a charge unit qo = eocto-E-o- 

For concreteness, we list here the values of these quan- 
tities for N2 at normal pressure, used in the simulations 



n 







2.3 /im 



1-12 



t w 3 • 10" " s , 
E ~ 200 kV/cm 



Qo 

Me 



: 7.56 • 10 7 cm/s 
4.7 • 10 14 e/cm 3 
i 380 cm 2 /Vs . 



(2. 



We now introduce dimcnsionless quantities by defining 



r = R a , t = t /t , 

q= («+ - n e ) e/q a , a = n e e/q 

J = ~je e/(q v ) , E = £/E . 



(2.9) 



Note that with our definition, j now plays the role of a di- 
mensionless charge current. If we furthermore introduce 
the dimensionless diffusion coefficient D as 

D = D e a /fi e E , (2.10) 

we obtain what we call the streamer equations 



d T a - V-j = a/(|E|), 

d T q + V-j = 0, 

q= V-E , 

j = a E + D Vcr 



(2.11) 
(2.12) 
(2.13) 
(2.14) 



where V denotes the gradient with respect to the dimen- 
sionless coordinate r, and where the "ionization function" 



/(|E|) = |E| a(|E|)/a 



(2.15) 



i s as sumed to vanish at zero field. Townsend's expression 
~q) yields: 



/ T (|E|) = |E| exp(-l/|E|) 



(2.16) 



In general, we will treat an ionization function with the 
properties J37| 



/(0) = = /'(0) , and /'(|E|)>0 for all |E| . 

(2.17) 



The dimensionless equations (2.11)-(2.14) now depend on 
only one internal parameter, the dimensionless diffusion 
coefficient D. For the values used in p|,p| for N 2 under 
normal conditions, D w 0.1, while according to the data 
given by Raizer S, for Ne and Ar, D w 0.3. We believe 
that typical values are gene rally in the range 0.1-0.3, since 
in the approximation (2/7), ao/Ea s=s I/e and since the 
ratio D e j '/j, e appears to be commonly of the order of volts 
for large fields, while / is typically of the order of several 
electron volts. 

We are now able, solely on the basis of the dimensional 
analysis above, to make a first semi-quantitative predic- 
tion about streamers. We will in practice be interested 
in external fields E + = 0(1) (for E + <C 1 and a^ 1 on 
the order of micrometers, the electron avalanche process 
becomes much too ineffective for streamer fronts to de- 
velop at reasonably small distances; also our scale separa- 
tion approach discussed in the introduction might break 
down). We can therefore expect that, for D values < 1, 
as is the case for N2 , front widths will be of order a^ , 
and that in addition the reduced electron density o~ far 
behind the front on the inner scale will be of order unity 
as well. This leads one to expect electron densities in the 
streamer body on the order of 10 14 cm~ 3 , in agreement 
with numerical findings. 



III. HOMOGENEOUS SOLUTIONS AND THE 
CONCEPT OF FRONTS 

A. Homogeneous States and their Stability 

The first task, when studying in general the propaga- 
tion of a front, is to identify the nature and stability of the 
states which the front connects. We expect the invading 
state, here the ionized one created by the front, to be sta- 
ble p8| , while the invaded state can in general either be 
metastable or unstable. Physically, we of course expect 
the non-ionized state to be unstable in a non-vanishing 
field in the present model. (In an attaching gas forming 
also negative ions, it is conceivable that the non-ionized 
state is metastable for not too strong fields.) 

Equations ( |2.11 )-(2.14) immediately yield, that sta- 
tionary homogeneous states simply are solutions of 



0-/CIE1) = 



(3.1) 



So, these stationary states decompose into two families: 
(i) Non-ionized states, with a = 0, E arbitrary: Since 
the density of free electrons vanishes, no ionization can 
occur, whatever the value of E is. If also the density of 
ions vanishes, V • E = 0. Since these states correspond 
to the physical situation far ahead of the front, we label 
them (+). Moreover, since we will need in particular the 



case in which the field ahead of the front is constant, we 
take E + =const. 

(ii) Completel y sc reened states, labeled (— ), with E = 
0, <7~ arbitrary J39[: Whatever the electron density, for 
E = impact ionization does not occur and thermal en- 
ergy is much too small to permit ionization. 

Since the steady states we consider as well as the 
equations of motion are translation invariant in space 
and time, the eigenstates of the linear perturbations are 
Fourier modes of the form 



Sa{r,t) 

<5E(r,i) 



Ei 



exp(ik • r + lot) . 



(3.2) 



We first investigate the linear stability of the non-ionized 
state er + = 0. Upon linearizing the equations about the 
zeroth order values (a + — 0,E + ), we find two branches 
of modes: 

(a) The first, trivial branch is a zero mode (u> = 0), 
with <7i = 0, expressing that the electron density re- 
mains zero. This zero mode accounts for the degeneracy 
of the non-ionized states, i.e., for the fact that there ex- 
ists a (+) stationary state for each value of E + . (For 
E + 7^ const., these zero modes express the degeneracy of 
all steady states with q + = V • E + for any ion density q + 
as long as the electron density a + vanishes.) 

(b) The second branch of perturbations is associated 
with fluctuations carrying a finite electron charge; its dis- 
persion relation is 



UJ 



+ = ik ■ E^ 



/(|E H 



Dk 1 



(3.3) 



withiw+k-Ei = (f(\E + \)-Lu+)a 1 . The first term on the 
r.h.s. of (3.3) simply expresses the fact that the electrons 
drift, to first order, in the electric field E + with velocity 
(— E + ). The real part diuj + , the sign of which determines 
whether fluctuations decay or are amplified, contains a 
destabilizing term, expressing that any small electron 
density fluctuation is amplified at rate /, and a stabilizing 
term, due to the diffusive spreading of electron charges. 
For k 2 < /(|E+|)/L>, dluj + > 0: non- ionized states are 
unstable against long-wavelength perturbations. 



We note, that the single Fourier eigenmodes (3.2) vio- 
late individually the physi cal constraint that a be posi- 
tive everywhere. But Eq. ( [3.3] ) also determines the time 
evolution of physically allowed fluctuations (wavepack- 
ets) that are superpositions of these eigenmodes. For ex- 
ample, one easily deduces from it Lozanski's expression 
J40f for the time-evolution of a Gaussian-shaped small 
electron density with arbitrary constants c±, c 2 > 0, 

( 5a(r,r)-c 1 e/(l E+ l) r (c 2 + ADt)- 3 ' 2 

e -(r + E+r) 2 /( C2 +4D7) ] (34) 

as long as linearization around the non-ionized state 
holds. As expected, the center of the spreading packet 
drifts with velocity — E + , while the total number of 
electrons it contains is amplified at rate / and the 
wave-packet stays Gaussian, with time-dependent width 



ci + ADt. Such ionization modes derived by linearizing 
around the non-ionized state are known as electron or 
ionization avalanches in the gas discharge literature. 

We now perform the same linear stability analysis for 
the completely screened states (a~ - const., E~ = 0). 
The fact that /'(0) = from Eq. ( |2.17| ) assures that the 
linear perturbations are not affected by ionization; the 
dynamics thus evolves with conserved particle densities. 

Again, due to the existence of a continuous family of 
screened stationary states, parametrized by a~ , the spec- 
trum contains a branch of w = modes. For the nontriv- 
ial branch, the dispersion relation is given by 



-Dk 2 



(3.5) 



while the eigendirection of such a perturbation is given 
by 



ci 



ik • Ei = <j\ + qi = 



(3.6) 



Since (a± + q\) i s th e dimcnsionless ion density of the 
linear mode, Eq. ( |3.6[ ) simply expresses the fact that ions 
are completely immobile in our model. 

Equation ( |3.5| ) expresses the fact that the completely 
screened (— ) states are stable, the decay of perturbations 
being due to the added stabilizing effects of overdamped 
plasmons (— a~~) and electron diffusion. The k — > limit 
of the plasmon mode leads to dielectric screening El) . 



B. The Mechanism of Front Creation 

Let us now investigate the dynamical evolution of an 
initial state in which the electron and ion densities vanish 
everywhere except in a small localized region. An exam- 
ple of such localized initial conditions is an initially Gaus- 
sian electron density, as in the simulations [g|js[| — under 
what circumstances initial conditions arc sufficiently lo- 
calized will become clear later. As long as the electron 
and ion densities are small enough, we can neglect in 
linear approximation the changes in the field as we did 
above when linearizing about the non-ionized state. As a 
result, both densities will grow due to impact ionization. 
If this were the only mechanism, the space charge would 
remain unchanged and the ionization would continue in- 
definitely. However, the electrons are mobile, and at the 
same time they start to drift in the direction opposite 
to the electric field E. If we neglect for the moment the 
diffusion, this drift has two effects: First of all, the elec- 
trons start to drift in the direction of the anode. Impact 
ionization then starts in previously non-ionized regions 
as well, so the ionized region expands towards the anode. 
Secondly, as the electrons drift while the ions stay put 
(on the fast time scale) , a charge separation occurs which 
tends to suppress the field strength in the ionized region. 
When the size of the initial perturbation and/or the time 
during which the avalanche has built up are large enough, 
the screening of the field becomes almost complete in the 



ionized region so that ionization stops there. The behav- 
ior in this region can be described by linearizing around 
the screened state as done above. After an electrically 
screened body of the ionized region has developed, the 
initial ionization avalanche is said to have developed into 
a streamer. Thus streamer fronts are strongly nonlinear 
and determined by two competing mechanisms, which 
dynamically balance each other: the ionization process 
which is strongest at the leading edge and the screening of 
the field due to the free charges which increases towards 
the rear end of the front. This balance also explains 
our finding that the ionization length and the screening 
length in the plasma behind the NSF are of the same or- 
der of magnitude for field values that are not too small. 
Technically speaking, the challenge in constructing the 
full front is to connect the two regimes linearized about 
the homogeneous states in an appropriate way through 
the nonlinear regime of the front. 

In the above discussion, we have neglected electron dif- 
fusion. In this case the NSF propagates towards the 
anode with a velocity that is at least the drift velocity 
of the electrons in the local electric field. The PSF, in 
contrast, is moving in the direction opposite to the drift 
of the only mobile species, the electrons. Its space charge 
is formed by the ions staying put, while the electrons are 
drawn into the ionized body. Propagation of a PSF is 
therefore only possible if the electron diffusion current 
overcompensates the drift current. This in turn implies 
that if the diffusion coefficient D is small, electron density 
gradients must be extremely steep. From this discussion 
it already becomes evident — and we will derive this be- 
low — that for an NSF, diffusion is a small correction 
for D -C 1, since drift and diffusion currents are acting 
in parallel directions. In PSF, however, diffusion has to 
overcome the drift, and as a result in this case the limit of 
vanishing diffusion is very singular. We will see in Section 
IV that this manifests itself through the emergence of a 
new inner length scale D/ao = D e /(ii e Eo), the diffusion 
length associated with the electron drift velocity. 

Of course, a charged front only screens the normal com- 
ponent of the electric field. This is why electric screening 
is efficient in the head of the streamer, while the field 
penetrates in the body of a single streamer in the sim- 
ulations f§||. Our planar front analysis thus serves as 
a first approximation for the mechanisms in the moving 
tip of the streamer finger. 



C. The one-dimensional Streamer Equations 

Let us now restrict our analysis to the case of plane 
fronts perpendicular to a constant electric field. Of 
course, in practice planar streamer fronts will be unsta- 
ble to deformations along the front (very much like in the 
Mullins-Sekerka instability in crystal growth PQ] ), but as 
explained in the introduction, the planar front analysis is 
a first step towards understanding the dynamics on the 



inner length scale a^ 1 and time scale to. As such, it is 
the first basic ingredient for deriving an effective interface 
model on scales ~^> a^ . 

We choose the x axis as parallel to the field and perpen- 
dicular to the planar front so that E = Eit and V = x d x . 
From the point of view of matched asymptotic expan- 
sions, the electric field in the non-ionized region before 
the front will vary adiabatically slowly on the "inner" 
time scale r of the front, the timescale on which the front 
propagates over a distance comparable to its width, be- 
cause the length scales of the outer problem determining 
the changes of E are assumed to be much larger than 
the inner scale a^ . For our study of the inner problem, 
we thus take the asymptotic field value E + in the union- 
ized region constant in time. Furthermore, we will use 
the convention that the unionized initial state into which 
the front propagates is at the right towards large positive 
values of x, so that there 

cr^cr+ = 0, q^q+ =0 , E -v E+ , d T E+ = , 

for x — > +00 , (3-7) 

which motivates now the use of the superscript +. We 
emphasize again that "x — > +00" should be interpreted 
on the length scale a^ of the inner problem in the sense 
of matched asymptotic expansions |13,nJ]. Far behind 
the front, i.e., for x — ► —00, the discussion of Section II 
leads us to expect a homogeneous stable state 

er^er-^0, q -> q~ = , E -> , 

for x — > —00 . (3.8) 



Which value o~~ will be dynamically selected and what 
the corresponding front velocity and profile are, for a 
given fixed value of the electric field E + before the front, 
is the selection problem, w e aim to solve. 

The boundary condition (3.7) allows an important sim- 
plification o f the equations in one dimension: If we insert 
( [2.13D into Q2.12D , we obtain 



d T E + j = h(r) 



(3.9) 



where h(r) is an arbitrary function of time which is con- 
stant in space. In view of the boundary condition (3.7), 
h(r) vanishes at x — > 00 and t hus ev erywhere. For planar 
fronts, the model Eqs. ( 2.11 )-( f04 ) then reduce to 



d T a = d x (a E) + D d x 2 a 
d T E = - o E - Dd x a , 



<rf{\E\) 



(3.10) 
(3.11) 



with space charge and electric current given by 

q = d x E and j = a E + D d x a . (3.12) 

We will refer to this set of equations as to the one- 
dimensional streamer equations. They are the basic equa- 
tions of this paper, on which the rest of our analysis will 
be based. 



Eq. ( J3.ll ) implies that the field decays behind the 
front, if no strong density gradients act against it. As 
we shall see later when we will discuss our simulation re- 
sults in Section VI, such strong density gradients often 
occur during the transient regime before a PSF emerges. 
Once, however, a front has approached an approximately 
uniformly translating state, the electron density o~~ be- 
hind the front is almost homogeneous and the field be- 
hind the front then decays to zero on a time scale 1/er - 
according to ( 3.11 ). Note, that the local decay of the 
field for any nonzero electron density is due to electrody- 
namics of conserved quantities that continues also after 
the impact ionization has been suppressed. 

We finally note that in the limit where the diffusion 
is small (D <g; 1), it is easy to identify the crossover 
time from the linear avalanche regime to that of streamer 
propagation in the case that the initial electron density 
is small and nonzero only in a very narrow localized re- 
gion. As explained in the beginning of this section, in 
the avalanche regime we can neglect the changes in the 
background field E + due to the build-up of the charges. 
The evolution of the electr on de nsity is then described 
by the linearized version of (3.10), a linear equation with 
drift, diffusion and growth. Hence, if the initial elec- 
tron density i s, e .g., Gaussian, the electron density will 
according to (3.4) remain a Gaussian profile, whose max- 
imum drifts with a velocity |E| in the direction opposite 
to the field and whose amplitude grows exponentially 
as exp(/(S + )r) . In other words, if the total initial 
electron charge is N e (0) = j dx a(x, 0), then the to- 
tal number of electrons in this avalanche regime grows 
as N e (r) = A e (0)exp(/(i? + )r). Likewise, the total ion 
charge grows exponentially, but if both the diffusion con- 
stant and the width and amplitude of the initial per- 
turbation are small, the electron drift will separate the 
negative electron charge and the positive ion charge al- 
most completely. The crossover to the nonlinear streamer 
regime will therefore occur when the total charge in the 
positively and negatively charged regions is big enough 
that screening of the field becomes appreciable, i.e., at a 
time r r when 



N e ( Tc )^\E+\ 



f(E+) 



ln[|£; + |/iV e (0)] 



(3.13) 



IV. UNIFORMLY TRANSLATING FRONT 
SOLUTIONS 

Above we already have introduced the idea, that fronts 
asymptotically approach some shape, which is indepen- 
dent of the initial conditions. This is based on our ex- 
perience pij-Pql with other examples of front propaga- 
tion into unstable states that the front will acquire some 
asymptotic shape and velocity in the long time limit, 
which will be the same ("universal") for a large class 



of "sufficiently localized" initial conditions that comprise 
most physically relevant initial states. This property is 
often referred to as the front selection problem. Our sub- 
sequent analysis will therefore follow the usual strategy 
in examples of this type: We will first show in this sec- 
tion that there generally is a one-parameter family of 
front solutions. In Section V we then summarize our 
present understanding of the front selection problem, and 
on the basis of this predict the properties of the selected 
streamer front. The numerical simulations that confirm 
our predictions are presented in Section VI. 

Uniformly translating fronts with velocity v are sta- 
tionary in a coordinate system moving with velocity v. 
If we denote this comoving coordinat e by £ = x — vt, the 
partial differential equations (pde's) ( 3.11 ) and ( 3.12 ) in 
this coordinate system become 



% a + <% (a E) + D d c 2 a + a f(\E\) , 



d T EL = vdtE-aE-D dt 



(4.1) 



A front translating uniformly with velocity v in the fixed 
frame x is stationary in this comoving frame, d T <rL = 
= d T E\c. As a result, the corresponding front pro- 
files are solutions of the ordinary differential equations 
(ode's). (We continue to use partial differential signs d$ 
even though the uniformly translating solutions are func- 
tions of the variable £ only.) 



Dd^c 



{v + E)d i a + ad^E + a f(\E\) 
D <%cr - v d^E + a E 



= , 

= o, 

(4.2) 



These equations are analyzed below. Both for D = and 
for D 7^ 0, they admit solutions for a range of values of 
the velocity, so we are indeed faced with the question of 
front selection. 

It is important to realize that not all the exact uni- 
formly translating front solutions of these ode's corre- 
spond to physically relevant solutions. In particular, 
any physical electron density a needs to be non-negative 
(c > 0), but as we shall see the set (4.2) admits PSF so- 



lutions where a goes negative. We expect these solutions 
to be unstable (in accord with the "nonlinear marginal 
stability" scenario Q), and also not te be approachable 
from an initial condition with a > 0. Hence they are nei- 
ther dynamically nor physically relevant. Furthermore, 
note that in our model the ion density qi (= p + a) can 



only increase due to impact ionization [Eqs. (2.11) and 
( 2.12 ) imply d T qi = o-f(E) > 0]. With our convention 
that the non- ionized state is on the right, this implies 
that uniformly receding front solutions with v < are 
unphysical. We will therefore call a uniformly translat- 
ing front solution physical if 



V > and <t(£) > for all £ 



(4.3) 



A. D = Front Solutions 

In contrast to the case D ^ 0, where we can de- 
rive properties of uniformly translating fronts only either 
qualitatively by discussing flows in phase spa ce or quan- 
titatively by numerical integration, Eqs. (4.2) for D = 
can be integrated explicitly. Doing so, we derive a sim- 
ple explicit expression for the electron density o~ behind 
the front in terms of the field E + before the front; this 
analysis generalizes an earlier result of D'yakonov and 
Kachorovskii p9fl , and explicitly illustrates the existence 
of a family of uniformly translating solutions. For NSF, 
these results extend smoothly to the case D ^ 0: The 
electron density a~(E + ) derived for D — will turn out 
to be a good approximation for D < 1, and the small 
overshoot of a at the rear end of the front visible in the 
three-dimensional simulations in Fig. 1(c), is also recov- 
ered for D = 0. For PSF, on the other hand, we will 
see that D acts as a singular perturbation, so that the 
class of D = PS'f-solutions that we derive here is not 
relevant for the PSF selection problem for D < 1. 

The ode's describing uniformly translating fronts for 
vanishing diffusion are found by putting D = in ( |4.2] ). 
These equations then become 



d ( 



{v + E)a 



*f(\E\), 



v d^ \n\E\ = a 



(4.4) 
(4.5) 



Upon in sert ion of the l.h.s. of Eq. (4.£) for a in the r.h.s. 
of Eq. (4.4), this equation can then be expressed as a 
complete derivative by writing 



ds 



(v + E) a + v 



dx 



/(M) 



(4.6) 



For physi cal fronts with v > and a > [see (4.3)], we 
see from (4.5), that E is a monotonic function of £, 



sgn (%£(£) J = sgn q(£) = sgn E(£) = sgn E + 

for all £ . (4.7) 

This allows us t o use E as a coordinate instead of £. Ac- 
cording to (|3.7[ ), before the front at £ — > oo the ele ctron 
dens ity vanishes, so a + = o-[E + ] = 0. Eqs. (4.6) and 
(4.7) together then determine a as a function of E as 



a[E] = 



E 



Pe+ [E] 



(4.8) 



with the function 

\E+ 



Pe+[E] 



dxlML = PlE+] {\E\] (>0). (4.9) 

\E\ X 



The function p E + [E] is nothing but the ion density, as can 
be deduced by inserting q = d^E into (4.5) and equat- 
ing the charge density q with p — a. The ion density 



p for D = turns out to be a function of E and E + 
only, and to be independent of the particular front shape 
parametrized by v. 

The fields a, p and E as a function of £ can be found 
by solving the implicit equation for E = E(£) 



d(. In |£| 



piWI^D 



E 



(4.10) 



which can be derived from Eqs. (4.5) and (4.9). 

Eq. (4.8) immediately shows that physically allowable 
solutions with a > and v > must have v + E > 
for all field values. Because of the monotonicity of E as 
a function of £, this is automatically satisfied for F5F 
with E + > 0, but for NSF this implies in particular that 
v+E + > 0; together with v > we thus have for physical 
fronts 



v > max 



0, -E+ 



(4.11) 



In physical terms, the condition v > —E + expresses, that 
the velocity of uniformly translating fronts must be at 
least the drift velocity — E + of free electrons in the lead- 
ing edge of the fro nt, where the field is strongest. (Re- 
member, that ( |4.7[ ) implies, that the field is monotonic 
in space.) 

the 



For all values of v o beyin g the inequality (4.11 
solutions of (18) and ( |f.lC| ) are proper, physical 



y al- 
lowable solutions for fixed E + ; within the context of the 
present model, this illustrates a general feature of front 
propagation into unstable states, namely that there exists 
a family of front solutions parametrized by the velocity 
l2t 



In Fig. 2(a), we plot the solution (4.8) for a as a func- 
tion of E for the fixed value of the velocity v = 2 in the 
case that the impact ionization fun ction f(E) is given by 
the Townsend expression Jt{E) of ( 2.16 ) as in the numer- 
ical simulations MM. Note that in this representation, 
the state behind the front at £ = — oo corresponds to a 
point on the a axis, and that the front solution <r(£), E(£) 
is represented in this diagram by the flow along one of the 
trajectories towards either the positive S-axis for PSF 
or the negative E'-axis for NSF for £ — * oo. Note further- 
more that a overshoots the value o~ (= er(£ — v — oo)) in 
the case of NSF. This property as well as the monotonic- 
ity of cr[E] and accordingly of o~(£ ) for positive fronts, 
follows immediately from Eq. ( |4.8| ). For NSF, it can 
also be observed in the three-dimensional simulations of 
Vitcllo et al. M, shown in Fig. 1(c). 

The smallest E + for which a front solution with v = 2 
is shown in Fig. 2(a), is E + = —1.999. For this value of 
E + , <r[E] continues to increase till E rs E + and then sud- 
denly decays to zero. A short analytical investigation of 
(18) shows, that this behavior develops into a discontinu- 
ity of a[E] at the point E = E + for v = -E + . a[E] then 
increases monotonically up to f(E + ) for E J, E + and 
then jumps to zero discontinuously at E + . This shock- 
like behavior stays unchanged under a parameter change 



to cr(£). It is further discussed and motivated in Section 
V. 

An immediate consequence of Eqs. (4.S) and ( f4.9| ) for 
the electron and ion density is that the value a~ behind 
the front (where E — > 0) is a simple function of the value 
E + of the field ahead of the streamer profile: 



o-(£H 



P\B+\[0] 



x 



(4.12) 



The virtue of this expression for the electron densit y g~ 
far behind the front as well as of the expression (4.S) 
for the ion density p throughout the whole front, is that 
it is independent of the velocity v, hence independent 
of whichever front profile is selected, provided that the 
limit D — > is smooth. We shall see later that this 
D = result remains relatively accurate for NSF fronts 
with D < 1, and compare it to the results of the simu- 
lations [Bpf in Section VI. For PSF, on the other hand, 
the above result will turn out to be less relevant due to 
the non-perturbative nature of the limit D — ► in this 
case. 

For the Townsend function f T (\x\) [Eq. ( |2.16| )1 the 
function a~(E + ) can be expressed as 



a-(E+) T = \E+\ E 2 [\E+\-^ 

= mie+D-eJie+I ■■') . (4.1..-;, 



where E n (z) is the exponential integral J43 ], 

We finally note that the second form of ( 4.13 ) shows 
that a" approaches /t for large fields, since /t S> E\ 
for E + S> 1. For E + of order unity, a and fx are 
still of the same orde r, a nd this shows (for small D) 



that the growth rate ( |3.3| ) of long wavelength unstable 
modes in the unionized state is comparable to the damp- 
ing rate ( |3.5| ) of stable modes in the plasma behind a 
NSF. For small fields, the strict bounds on E% |l| show 
that a~ ~ E + Jt{E + ), so that the approximate equiva- 
lence of these two time scales does not hold for E + <C 1 , 
but in the small field range our starting model is not very 
realistic anyway, because of the neglect of stabilizing re- 
combination terms. 



B. D / Front Solutions 

For D ^ 0, we can not obtain the uniformly translat- 
ing solutions analytically. Moreover, perturbation theory 
around the D = case is not simply possibl e, a s D ap- 
pears in front of the highest derivative in Eq. (12), so the 
diffusion term acts as a singular perturbation. As a con- 
sequence, Eqs. (12) reduce to a set of two coupled first 
order ode 's for D = 0, while three are required for D =/= 0. 
However, we can still easily demonstrate the existence of 
a one-parameter family of uniformly translating front so- 
lutions for D =/= through standard counting arguments 
for ode's. Building on the results of such an analysis, 
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the solutions can then be constructed by integrating nu- 
merically in a stable direction, using so-called "shooting 
methods" Q. 

To perform the analysis, it is convenient to write the 
equations as a set of three coupled first order ode 's. There 
is some freedom for the choice of the third variable: The 
standard choice would be a' = d^a, but for the discussion 
of the singular limit as well as for numerical stability, the 
charge density q has turned out to be the most convenient 
choice. The ode's ([12]) then become 



d^ a = 



aE — vq 
D 



g/Ogj) vE-vq 

v D 



(4.14) 



Just as we thought of the profiles for D = as describ- 
ing flow in a two-dimensional (a, E) phase space, we can 
now think of Eqs. ( 4.14J ) as describing a flow in a three- 
dimensional (a, E, q) phase space. The velocity v just 
plays the role of a parameter in the flow equations, while 
£ again plays the role of a time-like variable — see the 
sketch in Fig. 2(b). 

The steady states of the full pde's discussed in Sec- 
tion III correspond to fixed points of the flow: the point s 
(<t,0,0) on the <7-axis are fixed points of the flow (4.14) 
and correspond to homogeneously ionized plasma states, 
while the E-axis is a line of fixed points (0, E, 0) each of 
which corresponds to a non-ionized state with a = a' = 
and B^O. 

A uniformly translating front solution now corresponds 
to the existence of a trajectory in this phase space that 
starts at "time" £ = — oo on the cr-axis and flows to the 
E-axis for £ — > oo. The multiplicity of such solutions 
(i.e., whether they exist as discrete sets, or, e.g., as a 
one- or two-parameter family) can be determined as fol- 
lows. If we linearize the flow near an arbitrary point 
(ct~,0,0) on the cr-axis by writing (a, 0,0) = (cr~, 0,0) 
+j4exp(— A~£), we find the eigenvalue equation 



A - 



A- 2 -A-^-^) = 
D D 



(4.15) 



The fact that there is a zero eigenmode is a consequence 
of the fact that the cr-axis is a line of fixed points. For 
the two nontrivia l eig envalues (which correspond to the 
linearized modes (3.5) about the ionized and screened re- 
gion by equating ik ■ x = — A - and lu~ = A~v) we have 



A± = 



v ± y/v 2 + 4Do- 
2D 



(4.16) 



The eigenvalue AZ is positive, and hence gives a decaying 
exponential; thus points along the corresponding eigendi- 
rection flow into the cr-axis as £ increases. The eigen- 
value A~, on the other hand, is negative and hence cor- 
responds to an unstable eigendirection, with flow away 
from the axis. This implies that at each point (u^ 1 0,0) 



on the cr-axis, there is, for fixed v, a unique eigendirec- 
tion (— AZ, 1, AZ)Ei along which the flow is away from 
the axis. This flow can be followed in two anti-parallel 
directions, determined by the sign of E\ . The one flowing 
towards positive values of E is the beginning of a PSF 
front profile, the one flowing towards negative E is the 
beginning of an NSF profile. From these eigendirections, 
one derives that for PSF with field perturbations E\ > 0, 
the electron density decreases close to the a axis, while 
for NSF it increases. Accordingly, before reaching cr = 
for £ — ► oo, a NSF profile has at least one maximum of 
a, while a negative one can be (and is) monotonic. This 
generalizes our result for D = 0, and is consistent with 
the findings of Vitello et al. || shown in Fig. 1(c). The 
physical origin of the maximum of a in the rear end of the 
NSF profile is the screening of the field: Due to the low 
ionization rate in an already fairly suppressed field the 
ion density has already almost acquired its final value, 
so the electron density has to overshoot its asymptotic 
value a~ so as to make d^E < 0. The screening behind a 
PSF happens by suppressing the electron density faster 
than the ion density for increasing £, and so there a is 
monotonic. 

Let us now investigate the stability of the flow near a 
point (0, E + , 0) on the E-axis. Upon linearizing the flow 
equations (4.14) and writing the £ dependence of the per- 
turbations in the form exp(— A + £), we find the eigenvalue 
equation 



A 4 



A+ -A H 



E+ 



D 



f(\E + \) \ 
D J 



(4.17) 



Again, there is a zero eigenvalue due to the fact that the 
whole E-axis is a line of fixed points. The two nontrivial 
eigenvalues are 



A£ = 



E+ ± y/(v + E+Y - 4Df(\E^ 
2D 



(4.18) 



These e igenv alues can be related to ( |3.3| ) in the same 
way as fl4.15| ) could be related to (|3~5|). For v + E + > 0, 
the real parts of these eigenvalues are always positive, so 
that both eigendirections are stable. In other words, for 
E + > —v, all points near the E-axis flow towards this 
axis — in slightly more technical terms: there is a two- 
dimensional stable manifold flowing into each of these 
points on the _B-axis. For E + < —v, the flow is away 
from the _E-axis, and fronts with v + E + < can not be 



constructed. This generalizes (4.11) to D ^= 0. 

The existence of a one-parameter family of fronts with 
velocity v > —E + can now simply be understood as fol- 
lows. As we saw before, there is one unique PSF and 
one unique NSF trajectory flowing out of each point on 
the cr-axis for fixed v and D. Since the flow defined by 



Eqs (4.14) is continuous, we can expect each trajectory 
to extend smoothly p5[ . Once the flow gets near the 
-E-axis, we know from the above analysis that the trajec- 
tory will be attracted completely to the axis, provided 
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v is large enough. Thus, for each a~ and v, there will 
exist two unique trajectories, i.e., a unique PSF solution 
and a unique NSF solution. Since each of these trajec- 
tories flows into a unique point on the E-axis, the flow 
equations implicitly define a unique relation of the form 
c~ = o~ (v, E + ) for each of the two types of fronts. For a 
given value of E + , we thus have a one-parameter family 
of front solutions, parametrized by v. 

There are two important properties of the front solu- 
tions associated with their asymptotic l arge £ behavior. 
First of all, we note that according to ( 4.18J ) the eigen- 
values A± are only real for 



v > v* = -E H 



2y/Df(\E+\) 



(4.19) 



This implies that the corresponding front profiles can cer- 
tainly not approach the asymptotic state a = ahead of 
the front in a monotonic way for v < v*: When the eigen- 
values are complex, the front profiles have an oscillatory 
tail of the form exp[-(9?A + )£] cos[(3A+)£). Clearly, this 
violates the physical condition that the electron density a 
should remain positive, so solutions with — E + < v < V* 
are physically excluded: v* denotes, in the present case, 
the smallest velocity of physically allowable uniformly 
translating front solutions. 

The identification of v* as a bound on the velocity 
of physically allowed front profiles depends only on the 
structure of the eigenvalues A + associated with the linear 
flow near unstable states. There is a second, nonlinear, 
way in which the range of physically allowed values of v 
can be bounded. To understand this, note that for any 
v > v* , the asymptotic decay of o~(t;) for £ — ► oo for a 
uniformly translating profile will be 



a(0 = A_ e~ A -^ 



.4. 



-Aft 



h.o.t. 



(4.20) 



with real coefficients A- and A + . Here, h.o.t. stands for 
higher powers of the two exponentials generated when ex- 
panding the equations to higher than linear order in the 
variables. Clearly, the smallest eigenvalue A* governs the 
asymptotic decay of the profile provided A- ^ 0. That 
A_ will gcnerically be nonzero for an arbitrary velocity 
v follows again from the counting argument above for 
the flow in phase space: Each PSF and NSF trajectory 
flowing out of a point on the tr-axis is unique, and hence 
there is no freedom to impose an additional condition 
A_ = close to the E-axis. Furthermore, the coeffi- 
cients A- and A + depend on the full global nonlinear 
behavior of the flow, and hence they depend implicitly 
on v. 

There might exist, however, particular velocities 



,.part 



> v* , for which 



A_( v part )=0 



(4.21) 



For discussing these we invoke again a continuity argu- 
ment for the front properties as a function of v. Wc 
expect a very slowly decaying, nearly homogeneous uni- 
formly translating front solution to have a non-negative 



density everywhere, and to have a very large velocity, 
since the velocity of a profile is essentially inversely pro- 
portional to its slope in the li mit of small slopes. (So 
i ndee d the roots A I given by ( ]4.16| ) and A~t given by 
(4.18) vanish in the limit that v becomes large.) So 
for large v we expect to find physical solutions. These 
are characterized by A- > in the leading edge of the 
front. Decreasing v continuously, we either reach v = v* 
smoothly with still A- > 0, or we reach the first particu- 
lar velocity, v par , where A- vanishes. In the latter case, 
we expect by continuity A_(y) < for v < v pa7 . This 
implies that then a approaches zero from below, i.e., that 
the front solution is unphysical. Below the next v^ , we 
expect the electron density to develop two zero's and so 
forth. The largest v part , if it exists, thus plays the role 
of the nonlinear front velocity v^ |24j , 



max{/ ort \A„{v part ) = 0} 



(4.22) 



for a given E + . (Note that if A+ < 0.5A+, the higher 



order terms in ( 4.20D of order exp(— 2Al£) are actually 
larger than the second term exp(— A+£). This does not 
change our argument, though, as the prefactor of this 
second order term will vanish if A- vanishes.) 



At the velocity v 



or at any v = v part , the flow 



in phase space approaches the .E-axis along the eigen- 
vector where the flow is most rapidly contracting. The 
trajectory corresponding to the nonlinear front solution 
is therefore more appropriately referred to as a strongly 
heteroclinic orbit, where heteroclinic indicates that it is a 
trajectory from one fixed point to another one. The exis- 
tence and properties of strongly heteroclinic orbits have 
recently been under active investigation pG| . 

Such a velocity v\ if it exists, bounds the continuum 
of velocities of physical uniformly translating solutions 
from below, and thus replaces the earlier bound v* de- 
rived from linearizing the equations in the leading edge 
of the front. 



C. Nonlinear Front Solutions for PSF 



For NSF, the bounding velocity v* given by ( 4.19| ) is 
alway s pos itive. Moreover, by integrating the flow equa- 
tions (4.14) numerically and searching for particular so- 
lutions for which, according to (4.21), A-(v part ) — 0, we 
have convinced ourselves that there are no such solutions 
for any D ^ and E + < 0. Hence, the smallest velocity 
of physical NSF solutions is always v* , for any value of 
the parameters. 

For PSF, on the other hand, the situation is very 
different, since v* < for (E+) 2 > 4Df(\E+\) - for 
the Townsend function (2.16|), this happens for D ' 



O.ZhE+e 1 '^ , hence for any E+ for D < 0.68. In par- 
ticular for PSF at small D the question therefore arises 
whether there are nonlinear front solutions defined by 
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( [4.2l| ) and ( f4.22| ) with i; 1 * > 0. The results of a numer- 
ical search for such solutions are shown in Fig. 3, as a 
function of D and E + . Below the full line in this di- 
agram, there exists indeed a nonlinear front ifi > v*, 
whereas above this line v* denotes the smallest velocity 
of physical front solutions. While these results have been 
obtained numerically, the existence of a single (unique) 
particular solution with A- (v part ) = in the limit D — > 
can be demonstrated analytically. Since a full discussion 
of these results will be given elsewhere |47J], we confine 
ourselves here to a brief outline of the arguments that 
also demonstrate the singular nature of these solutions 
for D -> 0. 

If we take the limit D — > with v fixed, assuming no 
nontrivial scaling of the va riable s cr, E and q and of the 
spatial coordinate £, Eqs. (4.14) can easily be shown to 
reduce to those studied in Section IV. A for D = 0. Hence, 
we can recover in this way the family of front solutions ob- 
tained there. Any particular solution, on the oth er ha nd, 
for which A-(v part ) — 0, decays according to ( 4.20| ) as 
cxp(— A^£) as £ — * oo. Since A+ ex D^ 1 for D — ► 0, such 
a particular front solution becomes extremely steep as 
D — > 0: its gradients diverge as 1/D and so that the dif- 
fusion term can still overcome the drift term as D — > 0. 
That the velocity of such a solution must also have a 
nontrivial sc aling in this limit can be seen from the third 
equation of (4.14), written in the form 



d s q 



f(\E\) 



D q 



(4.23) 



Any nontrivial scaling of this equation in the limit D — ► 
can only occur if the first term between brackets remains 
of the same order as the other two, which diverge as 1/D. 
This is only possible if v scales as D. In this limit, the 
third term can then be neglected, and since d^q has to 
change sign in the front region (as the charge density q 
vanishes as £ — > ±oo), there must be an intermediate 
value E < E+ of the field for which v = Df(\E\)/E. 

Now that we know the scaling of the spatial gradient 
of the velocity of such particular front profiles for D — » 0, 
one easily convinces oneself that the electron and charge 
density of these solutions must diverge as 1/D in this 
limit. To study the existence of such possible solutions, 
it is therefore convenient to introduce new variables and 
coordinates according to 

X = Dx . v — Dv , £ = _D£ , a — a/D , q = q/D , 

(4.24) 

with E and r unch anged. In these new variables, the 
flow equations ( 4.14 ) become 



d;a = — a E + D v q, 
a ? ~ E = q, 



(4.25) 



The limit D — > can now be taken simply by leaving 
out the term Dvq in the first and last equation. We will 
show elsewhere 47 that the resulting equations have one 
unique physical front solution thus fixing one particular 
value of the scaled velocity v\ and in view of the scaling 
( 4.24 ) and the scaling of the eigenvalues A±, this solu- 
tion must have A_(v\) = 0. This solution is therefore 
precisely the D — > limit of the nonlinear front solution 
with velocity v^ = v\D. Furth ermore, since the limit 
D — ► is smooth for Eqs. (4.25), this shows that there 
exists a nonlinear front solution with ip > for any E + 
and nonzero but small D. Due to the singular scaling 
(4.24), the corresponding front solutions are determined 



by ode 's that have a different structure from those studied 
for D = in Section IV. A, and therefore these nonlinear 
front solutions can not be obtained perturbatively from 
the latter class of solutions — of course, the latter class 
of solutions still exists for D ^ 0, in agreement with the 
counting arguments given earlier, but these now corre- 



spond to a singular limit of Eqs. (4.25)! The significance 
of these nonlinear front solutions lies in the fact that 
they will turn out to be the selected fronts that domi- 
nate the dynamics of PSF in the physically important 
range 0.1 < D < 0.3. 

The nonlinear front solution can be constructed numer- 



ically very easily by integrating Eqs. (4.25) using stan- 
dard numerical "double shooting" routines J44|. Fig. 4 
shows our numerical results for the smallest physical ve- 
locity, max(ti*, v*) in the case that the ionization function 
is given by the Townsend expression. The scaled veloci- 
ties v' /D and v*/D are plotted; in agreement with our 
arguments above the scaled velocity v' /D of the nonlin- 
ear front solution approaches a finite limit as D — > 0. 
Furthermore, the ratio v^ / D hardly varies with D in the 
physical range 0.1 < D < 0.3, and for small fields E + , 
the scaled velocity v^ jD becomes exponentially small, in 
agreement with the bound v' /D < E + cxp(— 1/E + ) tha t 
follows from the observations discussed after Eq. (4.23) 
above. 

We finally note that our numerical routines have not 
only allowed us to obtain the results show in Figs. 3 
and 4, but have also enabled us to verify numerically 
all the statements made above about the multiplicity of 
solutions, the parameter ranges for physical fronts, the 
monotony properties, the singular behavior of the small 
D PSF limit, and the persistence of the family of front 
solutions for D — > 0. 



V. SELECTION OF THE ASYMPTOTIC FRONT 

A. Front Propagation into Unstable States 

We have seen that the non-ionized state into which the 
streamer fronts propagate, is an unstable state, that the 
homogeneous weakly ionized plasma is a stable state, and 
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that there is a family of uniformly translating front so- 
lutions connecting the two. The existence of a family of 
front solutions is a generic feature of front propagation 
into unstable states. We, therefore, briefly recall what 
is known in the literature for analogous problems and 
then translate this experience to the streamer problem. 
The prototype equation for studies of this type of front 
propagation is 



dtu = d^u + g(u) 



(5.1) 



where g(u) is some nonlinear function which satisfies 

ff (0) =0, 5 (1) = 0, </(0)>0, 5 '(1)<0. (5.2) 

Note that these relations imply that the "state" u = 
is unstable, and that the "state" u = 1 is stable. The 
study of the propagation of fronts into the unstable state 
u = in this equation dates back to the early work of 
Kolmogorov el al. p^| and Fisher p3 in the context of 
population dynamics. Later Gel'fand [|48| studied a par- 
ticular example of a function g(u) motivated by combus- 
tion. The mathematical research on this equation culmi- 
nated in the work by Aronson and Weinberger [Bp, who 



rigorously solved the front propagation problem for (5.1). 
In particular, they proved that any initial perturbation, 
that is nonzero only in a finite part of space, approaches a 
unique uniformly translating front solution with velocity 
Vf in the long time limit. If g"(u) < for all u, Vf equals 
v* = 2y // g'(0) (derived from linearizing in the tip of the 
front), while for general g(u), Vf approaches either v* or 
some v^ > v* . We refer to the literature for a detailed 
discussion of this work [J3l|]32j l . 

The velocities v* an d ir of the ab ove p roblem directly 
correspond to our v* ( 4.19 ) and v^ ( |4.22] ), since they are 



also the smallest velocities, which still allow for uniformly 
translating fronts with u > everywhere. So if u is in- 
terpreted as a population density or a chemical concen- 
tration, the selected front for every sufficiently localized 
initial state is the slowest physical uniformly translating 
front. In other interpretations no physical constraints 
bind u to positive values. Nevertheless the selected ve- 
locity stays the same. In this case, one can prove that 
every front with smaller velocity is dynamically unsta- 
ble plf, i.e., that the selected front is marginally stable. 
The slowest physical or stable solution, which is selected, 
coincides with the steepest physical or stable one. 

In the last decade, it has been recognized that sev- 
eral aspects of the front selection pr oble m encountered 
for the nonlinear diffusion equation (5.1), seem to have 
more general validity. Certain scenarios, justified by 
heuristic arguments but lacking a detailed mathemat- 
ical proof, were formulated and numerically tested on 
more complicated pde's that were often of higher order 
in the spatial derivatives [[49|]2l| -|25[| . Some of the equa- 
tions studied lead to non-uniformly translating fronts 
that leave a nontrivial spatially periodic state behind 
p9|,pl 23J24|,|5(J. A particular scenario is the one distin- 
guishing between the so-called linear marginal stability 



regime where Vf = v* and the nonlinear marginal sta- 
bility regime where Vf — v^> |2l|-f24||. These names stem 
from the fact that in this formulation, the two regimes 
of front selection are related to the stability properties of 
the front solutions — in both cases, the selected front sep- 
arates stable front solutions from unstable ones. Applied 
to (5.1), this scenario just provides an intuitive expla- 
nation of all the well-known mathematical results. For 
plasma physicists, it is worth mentioning that dynam- 
ics in the linear marginal stability regime is related to 
that determined by the "pinch point analysis" which was 
developed in plasma physics in the late fifties |5UB2L[23| . 



B. Predictions for Streamer Fronts 

By extending the arguments in the appendix of p3| , 
one may show that in the st ream er case just like in the 
case of the above problem (5T), all physical solutions, 
i.e., all solutions with u > resp. a > everywhere, 
are stable. For a detailed discussion, w e ref er to f47|| . It 
can be argued J22J23], and proven for (pTj) pi] , that a 
sufficiently localized initial condition will approach the 
physical uniformly translating front, which is closest in 



"phase space", i.e., the steepest one. Both for (5T) and 
for the streamer equations, the steepest uniformly trans- 
lating physical front is uniquely defined. It is also the 
slowest one. 

We can immediately prove this when initially u(x,t = 
0) = for x > x c for streamer fronts with D = 0: In gen- 
eral, there is a front solution for every v > max[0, — E + ], 
but now the only way in which the electrons can enter 
the range x > x c is through electron drift with velocity 
—E + . Clearly, therefore, the asymptotic front speed of a 
NSF can only be — E + , while a PSF can not propagate 
at all. If the initial electron density, however, decays ex- 
ponentially, the local electron density grows by drift and 
ionization, and the front can move quicker than — E + > 
for a NSF. 

For D y^ 0, we will here only conjecture the analogous 
statements, and we will test them numerically in Section 
VI: 

1. Selected front velocity. If the initial conditions 
are sufficiently localized, the selected front is the 
slowest physically acceptable front solution, i.e., the 
slowest front profile for which a(£;) > for all £. 
In view of the discussion of Section IV, this means 
that the selected front velocity Vf is predicted to 
be 



Vf 



■& 



2y/Df(\E+\) 



(5.3) 



except when there exists a nonlinear front solution 
satisfying (4.22): In that case 



Vf = v 1 



(5.4) 
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Note that the result ( |5.3j ) (v* is the Linear Marginal 
Stability value in the terminology of [|22|,|23|) is an 
explicit expression for Vf in terms of parameters 
associated with the linear instability of the unsta- 
ble state only. On the other hand, the existence of 
a nonlinear front and the value of v^ (the Nonlin- 
ear Marginal Stability value) depends on the w hole 
nonlinear behavior of the flow equations ( |4.14 ) . 



2. Localized initial conditions. Initial conditions are 
sufficiently localized if their spatial decay is faster 
than the asymptotic decay associated with the 
smallest eigenvalue of the selected profile, i.e., if 



a(x, t 



0)<Ce~ A », 

a(i,r = 0)<Ce 4 -(" t ) 1 , 



(5.5) 
(5.6) 



for 



depending on whether the selected front is v* or 
v* . Here C is an arbitrary constant, and A~t(v*) 



(= A~\_{v*)) and Al(t;t) are given by (4.1 



3. Non-localized initia l co nditi ons . If an initial condi- 
tion does not obey (|5.5|) or (^) , faster front speeds 
are possible. In particular, if initially a(x, t — 
0) ~ cxp(-Ax), with A < A±(v*) or A < At(u t ), 
whichever regime applies, then the front speed is 
given by 



-E+ + DA - 



f(\E + \) 



(5.7) 



which is obtained by solving (4.17) for v in terms 
of A. 

We now combine the analytic and numeric findings 
from Section IV with the selection rules above to quan- 
titative predictions for asymptotic fronts, which evolve 
from sufficiently localized initial conditions, in the case 
that the impact ionization is given by the Townsend ex- 
pression (2.16): 



NSF: For NSF, we numerically have not found any non- 
linear fronts for any D and E + , so our simple yet 
powerful prediction is that for NSF Vf = v* with 
v* given by d5.3| ). In principle it is possible that for 
other ionization functions f(E) than the Townsend 
function (2.16), there can be nonlinear front solu- 
tions also in the NSF regime. In practice, we ex- 
pect, however, that this will not be the case for 
physically reasonable funct ions f(E), i.e., for func- 
tions consistent with (2.17). 



Once the predicted velocities are known, the value 
(j~ of the electron density behind the streamer head 
is obtained from the numerical integration of the 
flow equations. The results of these calculations 
are shown in Fig. 5(a). Since for NSF, the limit 



D — > is smooth, also a depends only wea kly on 
D for D < 1, so that the D = prediction fl4.13j ) 
is quite accurate for realistic values of the diffusion 
coefficient. 

At the predicted values of the selected front veloc- 
ity, the width of the front region can be obtained 
directly from our numerical solutions of the flow 
equations. We have somewhat arbitrarily defined 
the width w as the distance between the points 
where a is 90% and 10% of the value a~ . As Fig. 6 
shows for NSF fronts with D — 0.1, this front 
width is typically of order 3 for field values of or- 
der unity. This confirms again that in the small 
D limit the impact ionization length a^ sets the 
inner scale of streamer fronts. Furthermore, we 
find that our numerical data are well fitted by the 
expression w rs 6/A±(v*), which shows that the 
front width simply scales with the spatial decay 
rate A+(v*) = A_(v*) of the streamer profile in the 
leading edge. NSF fronts always have a maximum 
of the electron density within the front. 

PSF: As we saw in Section IV, for PSF with D < 0.9, 
there always is a nonlinear front solution with ve- 
locity v' > v* . The prediction is that in this range 
the selected front solution is the nonlinear front so- 
lution, i.e., Vf — v^ . Values of v^ as as function of 
D and for several values of E + were already given 
in Fig. 4. We also saw before that these nonlin- 
ear front solutions are singular in the limit D<1, 
where u 1 * « Dv*(D = 0) and a~ « a~{D = 0)/D. 
The resulting predictions for o~~ are shown in Fig. 
5(6). 

The fact that the dimensionless inner decay length 
of these nonlinear fronts scales as D, implies that 
the physical decay length of such solutions is 
D/ao = D e /(/j, e E ), i.e., is given by the electron 
diffusion length. However, since simultaneously the 
electron density o~ diverges as 1/D, the total front 
width w defined above still approaches a finite limit 
as D — ► in units of the ionization length a^ . 

We finally note that the front propagation problem 
posed by the one-dimensional streamer equations has a 
number of interesting differences and similarities with 
the Aronson Weinberger front propagation problem ( f5.l| ) . 
In particular, it can be hoped, that techniques of strict 
bounds developed for the time development of these 
fronts |nj as well as for the nonlinear front velocity v^ 
p6[ might be also applicable to planar streamer fronts. 



VI. NUMERICAL TESTS OF THE PREDICTIONS 

We have tested the predictions listed in Section V by 
numerically integrating the pde's (3.1C), ( 3.1 1| ) forward in 
time. Our computer program is a finite difference code 
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with a time integration which is based on a semi-implicit 
method. 

We have performed an extensive search through pa- 
rameter space, varying D between 0.02 and 3, and \E + \ 
between 0.3 and 10. All our numerical studies of the 
dynamics fully confirm our predictions for fronts, and 
therefore we only present a sample of our results that 
illustrate the important features. 

All the simulations of the initial value problem we 
present in the remaining figures, have initially a field 
E = — 1 constant in space. We keep the field constant 
in time in the non-ionized region. The simulations of 
Fig. 7-10 start with the same localized initial ioniza- 
tion seed, a Gaussian profile for the electron density, 
a(x,t = 0) = 0.01 exp — (x — xq) 2 . Fig. 7 shows a run 
for D = 1 and times t = - 130 in time steps Ai = 2. 
As can be seen, the small ionization seed near xo = 50 
initia lly grows while drifting to the right in accord with 
Eq. (3.4). At time t = O(20), the ionization is strong 
enough that field saturation sets in and two asymmet- 
rically propagating fronts emerge. The one propagating 
to the right develops into a uniformly propagating NSF 
with velocity v* = 2.21 K| and degree of ionization be- 
hind the front o~ = 0.130. The maximum value of a 
in the rear part of the front is a ma x — 0.150. At the 
same time, a structure develops on the left, which at time 
t = 130 not yet has reached a stationary form, and which 
eventually will develop into a PSF . (Note that propa- 
gation to the left into a negative field — E + corresponds 
to a PSF front moving to the right towards x — > oo in a 
field +E + ). How the PSF actually reaches a uniformly 
translating profile, is shown in Fig. 8, where the devel- 
opment for xq — 150 and otherwise identical initial and 
boundary conditions is followed in time steps of At = 10 
during the time t = - 500. An asymptotic velocity of 
v< = 0.22 and a degree of ionization <j~ = 0.43 is reached. 
Note the huge difference in the degree of ionization and 
in the front velocity already for the unrealistically large 
value of the diffusion constant D = 1. 

The predictions from Section V for the selected uni- 
formly translating fronts for D = 1 and E + = ±1 yield 
for the NSF v* = 2.213 and <j~ = 0.129, and for the 
PSF «t = 0.2199 and a~ = 0.432. They thus correctly 
predict the simulations of the initial value problem shown 
in Figs. 7 and 8 within the accuracy given. Note that for 
the velocity v' of the PSF and for the degrees of ion- 
ization a~ both behind the PSF and the NSF this fact 
also shows the relative accuracy of the two very differ- 
ent numerical methods used, while for the velocity v* of 
the NSF the numerical integration of the initial value 
problem exactly reproduces the analytic result. 

As D decreases, both the structures within the fronts 
and the asymmetry between NSF and PSF become 
more pronounced. We illustrate this in Figs. 9 and 10 
with the temporal development starting from the same 
initial perturbation as before, but now for D = 0.1, 
the value corresponding to the parameter values of the 
earlier three-dimensional simulations MM. The time 



ranges in each plot are chosen appropriately for seeing 
the NSF and the PSF evolve into a uniformly translat- 
ing state. Fig. 9 shows a perturbation (initially localized 
at xq — 50) evolving during time 4 = 0- 130 in time 
steps Ai = 2. Except for the smaller diffusion constant 
and the stretched x axis, the situation is thus identical 
with that of Fig. 7. The NSF on the right propagates 
with a somewhat smaller velocity v* = 1.39, leaving a 
slightly higher ionization a~ = 0.147 behind. The max- 
imum <r max — 0.199 is relatively higher, since diffusional 
smoothening of structures is less pronounced. On the 
time scales of Fig. 9, the left front does not propagate, 
but retracts into an apparently immobile structure. The 
electrons drift with the field into the ionized region, leav- 
ing a layer of screening ions behind. Thus the electrons 
and the field are almost separated such that ionization on 
this side almost cannot occur. Eventually few electrons 
will reach the nonzero field region by diffusion and slowly 
build up a higher ionization and ultimately a propagat- 
ing PSF. That a PSF actually emerges, is shown in Fig. 
10. Only times t = 4000 - 8000 in time steps of Ai = 100 
after the initial perturbation at t = and xq = 60 are 
shown. The front propagates with velocity v' = 0.0149 
leaving behind an ionization o~ = 6.32. The numer- 
ical values predicted in Section V are v* = 1.384 and 
a~ = 0.144 for NSF, and v ] = 0.0146 and a~ = 6.234 
for PSF. The remaining numerical discrepancy of max- 
imally 2% could be resolved by choosing a still smaller 
gridsize in Figs. 9 and 10. Comparison of the PSF for 
D = 1 and D = 0.1 indicates that the time it takes such 
a front to build up, rapidly increases with decreasing D, 
but we have not pursued the scaling of the transient time 
with D. 

We finally show in Fig. 11 the evolution of streamer 
fronts starting from non- localized initial conditions, i.e., 
not obeying the bounds ( |T5| ) or ( |5.6P for D = 0.1. We 
used an initial electron density profile a(x,t = 0) = 
0.01/(2 coshA(x - 200)) with A = 0.25 and an ini- 
tial field E = -1. At these values, for the NSF, 
A±(v*) = 1.918 and for the PSF, At(w t ) = 0.3766. 
In this case, the bounds (5.5) or ( |5.6| ) are indeed vio- 
lated for both fronts, and Eq. (5/7) predicts a PSF with 
velocity v = 0.497 > w f = 0.0146 and an NSF with ve- 
locity v = 2.497 > v* = 1.384. The simulations find the 
fronts propagating with velocities 0.50 and 2.50, respec- 
tively. The ionization behind the NSF is a~ = 0.149 
and behind the PSF o~ = 0.158, so that now both are 
comparable to each other and to a~(D = 0) = 0.1485 
found from Fig. 5(a). Note that the diffusion constant is 
identical with that of Figs. 9 and 10, the only difference 
being the extended initial perturbation. 

The simulations confirm that streamer front propaga- 
tion is indeed correctly described by the marginal stabil- 
ity scenario, which in the present case amounts to the 
statement that the slowest physical velocity is selected, 
whenever one starts from sufficiently loca lized initial con- 
ditions, just as for the simpler case (pTj). 
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VII. CONCLUSIONS AND OUTLOOK 

The analysis in this paper fully supports the validity 
of an effective interface approach suggested by the re- 
sults of the full three-dimensional simulations of Dhali 
and Williams and of Vitello et al. [p|9]. This emerges 
from our detailed study of the associated one-dimensional 
problem, which yields the following results: 

(a) After a very brief stage of transient exponential 
amplification of the initial ionized seed, the grow- 
ing streamer evolves into an electrically screened 
plasma body separating two narrow fronts which 
propagate into the non-ionized outer region. We 
show that these two fronts correspond, for all prac- 
tical purposes, to translating profiles which prop- 
agate independently. This entails that the separa- 
tion of spatial scales between an inner front and 
an outer one, set by the global geometry, is indeed 
justified. 

(b) This enables us to draw upon the existing knowl- 
edge about front propagation into unstable states 
and thus to provide definite predictions about: 

• the relationship Vf(E + ) between the velocity of a 
planar streamer front and the value of the electric 
field ahead of it, and 

• the value of the degree of ionization of the plasma 
created by the front, a~{E + ). 

These predictions, although only valid as such in 
the absence of front curvature, still compare very 
favorably with the numerical results of Ref. ||. 
The two values of a~(E + ) on the axis of Fig. 1(a) 
and 1(6) behind the curved fronts of the 3D sim- 
ulations M (with the convention that E + should 
be understood as the electric field value extrapo- 
lated from the external non-ionized region to the 
front position) are plotted in Fig. 5(a). Without 
adjustable parameters our one-dimensional predic- 
tions for a~(E + ) are well within a factor of 2 from 
the 3D simulations. Likewise, the velocity values 
for vj(E + ) even agree to about 20%. 

Moreover, our analysis shows that NSF and PSF 
propagate in this model and for realistic values of the re- 
duced diffusion coefficient D, in a very asymmetric man- 
ner: 

• NSF rapidly reach a regime of uniform propagation 

- typically, on the scale of several tens of time units, 
i.e., in less than 10 _10 s. Their velocity is slightly larger 
than the electron drift velocity in the field E + . 

• This is to be contrasted with the dynamics of PSF: 
For realistic D-values, of order 0.1-0.3, they approach 
uniform translation considerably more slowly than NSF 

— typically on the time scale of 10 _8 s. Moreover, their 
asymptotic velocity is also much smaller than vV SF '. It 
obeys the inequality vf SF < DE+ cxp(-l/£ ,+ ) g7J- Fi- 
nally, while the widths of PSF and NSF are comparable, 



the degree of ionization behind PSF is much larger (up 
to a hundred times for D = 0.1) than that behind NSF. 
These results answer the question of whether PSF do 
or do not propagate, while explaining why the simula- 
tions of Vitello et al. J9| could not yield a definite answer 
— most probably, because, although their total width is 
of order a$ , their true inner length scale, as defined by 
the steepness of the profile, was too small to be resolved 
by their grid size. (Note that the apparent symmetry be- 
tween PSF and NSF found in earlier simulations || is 
to be related to the fact that there propagation into a pre- 
ionized medium (with initial electron density of 10 8 /cm 3 ) 
is studied, and possibly also due to the use of a poorly 
resolving grid.) 

It was observed empirically by Dhali and Williams || 
that in the three-dimensional simulations, the dielectric 
relaxation time in the plasma behind the front was of 
the same order as the intrinsic time scale set by the front 
motion. Our analysis shows, that this was no accident: 
It is a manifestation of the fact that the balance of the 
growth mechanism (impact ionization) and the stabiliza- 
tion mechanism (screening) leads to a single time scale 
to = (aoA'e-E'o) -1 for a NSF and for the relaxation be- 
hind it for fields of order Eq. Since our dimensionlcss 
value of a~ is the inverse dielectric relaxation time, it is 
of order unity (or slightly smaller) for fields E + s=s — 1. 

Of course, the above results should only be considered 
as a first step towards a realistic treatment of streamer 
propagation. They will have to be developed and ex- 
tended along two different directions: 

(i) Predictions of patterns within the present model 
and comparison with the simulations: Within the frame 
of the present continuous and fully deterministic model, 
we here have only considered the restricted case of a 
one-dimensional geometry. This enabled us to demon- 
strate that the concept of effective interfaces does apply 
to streamers. This approach will now have to be extended 
to the description of curved fronts. As also discussed in 
p9| , one will then be equiped with a reduced formulation, 
valid on the outer scale, which will permit us to study 
real three-dimensional streamers as pattern-forming sys- 
tems, as was done, e.g. , for viscous fingers and dendritic 
solidification fronts |2(J • This should provide a direct ap- 
proach to the question of dielectric patterns, alternative 
to the phenomenological Z?L^4-inspired dielectric break- 
down models |B4J. 

(ii) Possibly, extensions of the model will be necessary 
to predict real experiments: We have based our analysis 
on the minimal model as defined in Section II. It contains 
several restricting simplifications. A first step in the im- 
provement of the model would be to include the field 
dependence of the transport coefficients D e and fi e . It is 
clear that this will not modify our qualitative analysis, 
as, e.g., the counting argument for the existence of front 
solutions in Section IV depends only on the lineariza- 
tion about the stable and unstable states. Moreover, the 
qualitative asymmetry between the NSF and PSF will 
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persist as these result from the asymmetry of the electron 
drift. Quantitatively, the value of v* , the selected value 
of NSF : will simply be given by ( p73j ) with the trans- 
port coefficient and ionization rate evaluated at the field 
value E + . The slow transient build-up and small speed 
of PSF could be affected quantitatively by ionic motion, 
but from this effect, we expect no major qualitative dif- 
ferences. 

Finally, it should be kept in mind that our continuum 
equations are only valid on length scales larger than the 
mean free path £ m f p . On the other hand, we find for the 
strongest field values appearing in the simulations (which 
are much larger than the values of the field across the gap, 
due to the enhancement near the streamer tip), that the 
front width decreases down to about 3o;q « 3€ m / p in 



the approximation (2.7). In such limits, nonlocality of 
the transport and ionization effects begin to play a role. 
In addition, under these conditions, a typical volume of 
size V? m t v contains only of the order of 1000 electrons 
for the parameter values (p^q) used in the simulations. 
Fluctuations are then likely to become non-negligible. In 
principle, treating these effects calls for a full kinetic de- 
scription. This is probably out of reach for the moment, 
but one might want to mimic the main features of these 
effects by introducing stochastic terms in the equations. 
These also could mimic photo-ionization somewhat be- 
fore the front due to photons released in the impact ion- 
ization events, or the natural homogeneous background 
ionization due to radioactivity and cosmic radiation. In- 
vestigation of their relevance for branching of dielectric 
breakdown patterns might help to understand the asym- 
metry between the macroscopic patterns of discharges 
propagating into a positive or a negative field p5[ . 

In conclusion, our analysis opens the way to a micro- 
scopically based interface approach to discharges that 
seems promising for building a coherent framework for 
the analysis of breakdown patterns of various degrees of 
complexity. 
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APPENDIX: DIFFERENCES AND 

SIMILARITIES BETWEEN COMBUSTION AND 

STREAMER FRONTS 



In the introduction, we draw on the similarity between 
the streamer problem and other problems like combus- 
tion, chemical waves, thermal plumes, phase field mod- 
els, etc., to motivate the development of an effective in- 
terface approach. Of these problems, streamer propa- 
gation is most closely analogous to combustion, in that 
the strong nonlinearity of the reaction rates (the combus- 
tion rate and the ionization rate) is an important factor 
in giving rise to front development in flames and stream- 
ers, respectively. There are important differences as well, 
however, and since several interface techniques were orig- 
inally developed in the context of combustion [0,ll2|,|l7f , 
we highlight some of the differences and similarities here: 

(a) In combustion the reaction rate depends strongly 
on the temperature, whose outer dynamics is governed 
by a diffusion equation of the form dtT = V 2 T, while 
for streamers the ionization rate depends strongly on the 
field |E|, with E the gradient of the potential $ that obeys 
the Laplace equation V 2< & « in the outer region where 
the total charge density vanishes. This field strength E 
varies strongly in the streamer front, since the increased 
screening resulting from the rising electron density sup- 
presses E — and hence the ionization rate — to zero. In 
combustion, on the other hand, the temperature hardly 
varies throughout the combustion zone. 

(b) Combustion fronts are essentially fronts progating 
into a metastable state, because the front has to supply 
the heat that increases the temperature and hence the 
reaction rate, while streamer front propagation is an ex- 
ample of front propagation into unstable states, where 
the reaction starts for any nonvanishing electron density. 

(c) In a flame front typically the temperature remains 
high enough that all the reactions proceed to saturation: 
all the combustable material burns. The temperature 
difference between the flame front and the background 
is then essentially determined by conservation (conver- 
sion) of energy. In typical streamer fronts, on the other 
hand, the field E is suppressed long before saturation ef- 
fects start to play a role, and hence the ionization level 
behind the front is set by the internal dynamics of the 
front rather than by conservation laws (i.e., the gas den- 
sity). 

(d) The electron drift — /i e E has no clear analogue in 
combustion. 

(e) Finally, the relevant asymptotic expansion for 
streamers is not quite like the "activation energy asymp- 
totics" of combustion (ljjl^l , since we consider here fields 
strengths that are comparable to the characteri stic field 
scale Eq of the ionization rate given in Eq. ( |2.6| ) be- 
fore the front, whereas in combustion activation energy 
asymptotics is often appropriate since the flame temper- 
ature remains much smaller than the chemical activa- 
tion energy. For streamers, an analysis like activation 
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energy asymptotics is appropriate in the limit of small 
fields |E| -C -Bo- Of course, in streamers the rapid vari- 
ation of the field E in the front region, and hence the 
rapid suppression of the ionization rate, looks, at first 
sight, similar to the suppression of the chemical reaction 
rate with decreasing temperature in flames. However, in 
flames this is due to the strongly nonlinear dependence 
of the reaction rate on temperature before the front, (so 
that a slight suppression of the temperature reduces the 
reaction rate dramatically), while in streamers in large 
external fields of order E this is due to the fact that the 
field itself is reduced significantly behind the streamer 
front, as a result of screening. 
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FIG. 2. (a) Uniformly translating fronts for D — and 
w = 2 shown as flows in the two-dimensional (E,a) phase 
space. Out of each point a~ on the a axis, there is a PSF 
flowing to the right and a NSF to the left. Both reach 
the same value |-E + | on the horizontal axis, which also is 
independent of v. Note, that NSF have a maximum of 
a within the front, while PSF have monotonic a. Note 
also, that no physical fronts (i.e., with a > everywhere) 
r each a value E + < — v — —2, in agreement with Eq. 
(4.11). (b) Sketch of a uniformly translating PSF and 
NSF for fl/Oasa flow in three-dimensional (E, a, q) 
phase space. The thick curves indicate the trajectories, 
while the thin ones show their projection into the a — 
and q — planes. For fixed v, there is at each point of 
the a axis still only one outgoing vector, which can be 
followed in two antiparallel directions. The E axis is fully 
attractive and will always be reached. 



FIG. 1. Results of the numerical simulations of the 



full three-dimensional streamer equations (2T)-(2.6) of 
Vitello et al., reprinted from Fig. 1 and 10 in M. (a) 
Negative streamer propagating downwards towards the 
anode. Electrodes are planar and located at z— and 0.5 
cm; the voltage between the electrodes is 25 kV, which in 
the absence of the streamer amounts to a constant elec- 
tric field \£\ — -Bo/4. The system continues sidewards 
sufficiently far to make the lateral boundaries irrelevant. 
The streamer is assumed to be cylindersymmetric. The 
dimensionless diffusion constant is D = 0.1. Each line 
indicates an increase of n e by a factor 10; densities of 
10 3 - 10 14 , 
tion: 1 cm 

ionization seed was placed near the upper electrode. (6) 
Shape at time 5.5 ns. (c) Logarithmic electron n e and 
total charge n 3 density along the symmetry axis of (&). 
Solid line: n e ; dot-dashed: \n 3 \ for n 3 > 0; dotted: \n s \ 
for n s < 0. Note the exponential increase of the densi- 
ties on fim scale within the front as well as maximum of 
both densities in the rear part of the front. Courtesy of 
P.A. Vitello. 



u 3 can be seen (initial background ioniza- 
3 ). Shape at time 4.75 ns after an initial 
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FIG. 3. Phase diagram for PSF as a function of D and 
E + . Above the solid line the lowest speed of physical 
front solutions is given by v* , below the line v' corre- 
sponds to the smallest speed of physical front solutions. 
Accordingly, the selected front speed is v* above the solid 
line (linear marginal stability regime), and v* below the 
solid line (nonlinear marginal stability regime). The dot- 
ted curve indicates v* — and is a lower bound for the 
cross-over to it behaviour of the selected fronts. 
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v*forE + = 




FIG. 4. v ] = v^/D (solid) and v* = v*/D (dashed lines) 
as a function of D for E + — 0.3 - 1.0 in steps of 0.1, 
and for E + — 1.0 - 2.0 in steps of 0.2. v' depends only 
weakly on D, i.e., the physical front velocity v' is ap- 
proximately proportional to D. At v^ (E + , D cr (E + )) = 
v*(E + ,D cr (E + )), the selected front crosses over from v' 
to v*; the v* solutions disappear. Plotting D cr (E + ) in 
the (E + ,D) plane yields the solid curve in the phase dia- 
gram of Fig. 3, while the zeros of v* determine the dotted 
curve in Fig. 3. 



FIG. 5. Electron density o~ behind the planar selected 
front as a function of the field E + before the front for sev- 
eral D\ dotted: v* fronts; solid: v* fronts, (a) NSF: For 
v* fronts, a~ depends but weakly on D. Results for D = 
0, 1, 3 are shown. Crosses: Extrapolation of a~(E + ) for 
D = 0.1 for the curved fronts of the 3D simulations M. 
(b) PSF results for D = 0.05, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 
2, 3. For v f fronts, b~ = const. + 0{D), i.e., a' (E + ) is 
approximately proportional to 1/D. 
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FIG. 6. Width w of the front profile (measured between 
points with 0.1 o~ and 0.9 a~) as a function of E + for 
the selected NSF fronts with D = 0.1. The dashed line 
is given by w = 6/Al (v*). 
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FIG. 7. Num erical i ntegr ation of the time evolution given 
by Eqs. (3.1C) and (3.11) for D = 1.0 in a constant back- 
ground field E — — 1 (numerical gridsize Ax = 0.1 and 
timestep At = 0.05, initial perturbation at xo — 50). Ini- 
tial condition at t — 0: small charge-neutral, ionized re- 
gion of Gaussian shape depicted by the lowest line. Each 
new line corresponds to a time step At = 2 and the upper 
line to i = 130. (a) The electron density <j(x, i) initially 
grows and then, after field screening in the middle sets 
in, develops into a NSF propagating to the right and a 
PSF propagating to the left, (b) The electric field E(x,t) 
stays E = — 1 in the non- ionized region and becomes dy- 
namically screened to zero in the ionized region. 




FIG. 8. Emergence of the uniformly translating PSF on 
the left in the system of Fig. 7. Conditions identical to 
Fig. 7 except for xo = 150 and different numerical grid- 
size (Ax = 0.05 and Ar = 0.01). <r(x,t) is shown in 
time range t — - 500 in time steps At = 10. Note 
the difference in the duration of the transient regimes, in 
the propagation velocities of PSF and NSF, and in the 
degrees of ionization behind these. 




FIG. 9. Identical with Fig. 7(a), except that here D = 
0.1. Time range also t — - 130 in steps of At = 2. The 
NSF has sharper contours and propagates slower than 
for D = 1, the PSF appears not to develop. 




FIG. 10. Emergence of the uniformly translating PSF 
on the left for D — 0.1. Initial conditions identical with 
Fig. 9. The time range t = 4000 - 8000 after an ini- 
tial perturbation at t — and xo = 60 is shown in time 
steps of At = 100. (Numerical gridsize Ax = 0.01 and 
At = 0.5.) 
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FIG. 11. A non-localized initial condition with A = 0.25 
as described in the text; otherwise, the situation is like 
in Fig. 9, and D = 0.1. 
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